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ABSTRACT 

An implicit approximate factorization (AF) algorithm is constructed which has the following 
characteristics. 

• In 2-D: The scheme is unconditionally stable, has a 3 x 3 stencil and at steady state has a 
fourth order spatial accuracy. The temporal evolution is time accurate either to 1st or 2nd 
order through choice of parameter. 

• In 3-D: The scheme has almost the same properties as in 2-D except that it is now only 
conditionally stable, with the stability condition (the CFL number) being dependent on the 
"cell aspect ratios,” Ay/Ax and Ax/Ax. The stencil is still compact and fourth order accuracy 
at steady state is maintained. 

Numerical experiments on a 2-D shock-reflection problem show the expected improvement over 
lower order schemes, not only in accuracy (measured by the Li error) but also in the dispersion. 

It is also shown how the same technique is immediately extendable to Runge-Kutta type schemes 
resulting in improved stability in addition to the enhanced accuracy. 1 
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under USAF Grant No. AFOSR-87-0218. 


1 INTRODUCTION 


It can be shown [1] that numerical approximations to the linearized Euler equations of gas dynamics 
give rise to dispersive errors which in the 2-D supersonic case depend on a similarity parameter 
k = (Ay/A x)\/M* 0 - 1 (under the assumption v <SC u everywhere, where u and v are the x 
and y components of the velocity vector). The difference between the dispersion relations of any 
numerical algorithm and that of the original partial differential equations can be plotted as curves 
in the Fourier plane with k as a parameter. 

In particular the results in [1] indicate that for central-difference schemes, the dispersive error 
are contributed mostly by the third power of the errors in the Fourier variables 6 and <j>. It is, 
therefore, natural to think of fourth order spatially accurate algorithms as having better dispersive 
properties. By utilizing the structure of the Euler equations one can obtain, on a cartesian grid, a 
fourth order approximation which instead of using a 5 x 5 stencil (and 5 x 5 x 5 in 3-D) relies on 
a compact support of 3 x 3 (and 3 x 3 x 3 in 3-D). The advantages of the combination of fourth 
order accuracy together with compact support are quite obvious in terms of total computer work 
and memory. 

In Section 2, we describe the construction of an approximate factorization (AF) central differ- 
ence scheme and examine its theoretical (linear) stability properties. In Section 3, we derive the 
corresponding 4-step Runge-Kutta scheme and show how the Jameson-Schmidt-Turkel algorithms 
[2] may be easily modified to that form which has, in addition to the higher accuracy, markedly en- 
hanced stability limits. In Section 4, we describe some numerical experiments using the AF-version. 
Section 5 summarizes our findings. 

2 DERIVATION OF THE APPROXIMATE-FACTORIZATION 
SCHEME 

2.1 The Two-Dimensional Case 

Consider a general hyperbolic conservation law in 2-D: 

+ /* + 9y= 0. (1) 

In the case of Euler equations, for example, the vectors u, / (u), 5 (u) are given by 
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where p, u, v, E, and p are respectively the density, velocity components in the x and y directions, 
the total energy per unit volume and the pressure. In addition there is the equation of state relating 
algebraically (in the case of gas-dynamics) the internal energy to the pressure and density. One 
may also write (1) in non-conservation form as 

u t -\-A u x -\-B Uy= 0 (3) 

where A and B are the Jacobians of / and 9 with respect to u . 
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Consider first a forward-time second order central finite difference approximation to (1) 
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where we have dropped the sup-arrow indicating a vector and used the conventional differencing 
notation 

Vj,k= u y,fc = *Ay, nAt). (5) 

We have a cartesian grid with nodes at x — j Ax, y = kAy and t = nAt. We shall now introduce 
the usual shift operations notation: 
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Eq. (4) then may be written as: 
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where A = At/ Ax and 2R = Ay/ Ax. If we take a Taylor series expansion about the grid point 
(x,y, t) = (j Ax, kAy, nAt), we can construct the modified equation corresponding to (7): 


At 

Ut + ~2 Uu + 



' Ax 2 ' 


Ay 2 

*— 

JX H ~~JXXZ + * * * 


9y + -3} ~9yyy H 


or 


u< + /* + 9y — 


At Ax 2 Ay 2 1 

T u “ + IT'*** + ~W 9vvv ' 


( 8 ) 


Thus if we want to approximate (1) to a higher than second order (and in particular to fourth 
order in space-see comment in Introduction concerning dispersive errors) we must modify (7) by 
subtracting out the terms on the right hand side of (8). At this point we realize that using a 
straightforward approach to approximating f xxx and g yyv will lead to bigger stencils. However, 
using the original partial differential equation, (1), we have 


fxxx — u txx 9 yxx 


(9) 


and similarly 


9\ tw — 
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( 10 ) 


Since f xxx and g yyy need be approximated only to second order (because of their coefficients, Ax 2 /3 
and Ay 2 /3) the required stencil for the terms on the right hand side is only 3x3. (This observation 
was previously made, in another context by Jones, et al. [3].) So, replacing f xxx and g yyy in (8) 
by -6%Dtu/Ax 2 At- 6 2 n y S y g/ Ax 2 Ay and —S 2 Dtu/ Ay 2 At — 6y 2 (t x 6 x f/Ay 2 Ax, respectively; and 
also replacing u« by D t u t /At — *• D t (-f x - 9y)/At -D t [(n x 6 x f / Ax) + {n y 8 y g / Ay)\ we obtain 
the following approximation to (1) which is spatially fourth order accurate and temporally second 
order accurate: 
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Using (3) and the definition of D t we rearrange the above into the delta form 
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where I is the identity matrix. This is an implicit unfactored scheme. Since we are interested in 
marching to steady state we approximate-factor the left hand side to obtain: 


(/ + + §^*a) (/ + + !*!'■•*») (“S 1 - <*) 

= -A (/+ i«;) (/+ !«]) ,?»} . ( 12 ) 

We introduced into (12) the parameter a. If a = 1 then we retain the temporal second order 
accuracy while a = 2 gives first order in time. Note that (12) involves the inversion of block- 
tridiagonal matrices. The right hand side represents the steady state operator to fourth-order. The 
whole scheme involves only a 3 stencil. 


2.2 The 3-D Case 

The starting point, corresponding to (1) is 

+ /* + g v + h z = 0. 

Following the same steps as in the 2-D case, we get the modified equation 
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where using (13), we have 
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Repeating all the steps leading to (12) we obtain its three dimensional analog: 

(a+^ + |aam.) ('+***? + § a |m») (' + ip J + f A i' ,A ) m) = 

= - a [*a (/ + + is;) n Kl + „,s, (/ + is; + is;) + „,s, (/ + is; + is;) *? J . 

( 18 ) 

In (18) the matrix C is the Jacobian d h /d u and the shift operators fi Z} S z are defined in a 
manner analogous to (6). Q is the cell aspect ratio in the z direction, Q = Az/Ax. 
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2.3 Stability Analysis for the 2-D Case 


We consider the linearized (i.e., frozen coefficients) version of (12). We carry out a Von-Neumann 
analysis in the usual manner by casting (12) in Fourier space. A typical Fourier component of u ” k 
is given by 
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where 


<f> = lAx, 9 = mAy 


( 20 ) 


with 


— tt < <f>,9 < ir and — oo < £, m < oo. 
The mapping of the various shift operations is as follows: 


( 21 ) 


S t ul k —* 2 isin$- = 2 »£, —* 2»sm| = 2try (22) 
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S l u ?,k-> -4«m 2 | = -4£ 2 , £ 2 u£,. -► -4sm 2 |= -4y 2 (24) 

and so, 

M**zti£* -► 2iCs/l - £ 5 , -» 2»r?y / l ~ »? 2 - (25) 

We also define the amplification matrix G by ti n+1 = Gti n . With these notations, the frozen 
coefficient version of (12) is mapped into the Fourier space as follows: 


(/ _ l A ? + tVAA^l - e 2 ) (/ - §2*1* + ia^y/l^Pj (G - 7) 

= -2, -A [a^/i - e (/ - In 2 ) + ^n\fi - n 2 (/ - §£ 2 )] • (26) 

In the general case the matrices A and B do not commute, thus rendering the analysis of (26) 
almost intractable. It is instructive, however, to consider the scalar case. Since the aspect ratio 
1R = Ay/ Ax is arbitrary, and since the original partial differential equation (1), in the scalar linear 
case, could be transformed to the wave equation 


u* + U| + up = 0 


with t = At, x = x and y = Ay/B, we can without loss of generality (in this special linear scalar 
case) rewrite (26) as: 
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Equation (27) can be rearranged to solve for the amplification factor G : 
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where 


(29) 


o. = \ — -£ 2 — -17 2 + -£V - — CT 2 A 2 ^T7^/l - ?\J\ ~ rj 2 
and b is the Fourier map of the steady-state operator, i.e.: 

4 = 2 [a^I - f>(l - ?„’) + - §£’)] ■ (30) 

We see immediately from (28) that for a = 1 (i.e., the case of second order temporal accuracy) we 
have a Crank-Nicholson type scheme with ||G|| == 1. For a > 1 (i.e., first order accuracy in time), 
we have ||G|| < 1. We have thus demonstrated the unconditional stability of (27) for all values of 
the cell aspect ratio. 

2,4 Stability for the 3-D Case 

The three dimensional analog to (27) is 

(1 - + ' a A £\/i - £ 2 )(i - |>? 2 + ~ ^K 1 ~ p 2 + ~ f2 )( G - x ) = 

W *(1 - §e 2 - |» 2 )' 
(31: 

where f is the Fourier dual variable defined analogously to £ and ij. 

The stability of three dimensional amplification factor G, as defined by (31), is difficult to 
analyze and we resorted to numerical evaluations of |G| 2 , using up to 8 x 10 6 Fourier modes. The 
numerical study of (31) was carried out on the Cray 2. We found that as in the case of forward 
Euler approximate factorization second order scheme [4], the amplification factor was conditionally 
stable. For example for JR = Q = 1, the stability limit is A < .43. These stability properties are 
obtained with the aid of artificial viscosity (AV) term which is of sixth order but still resides on 
the compact stencil. The AV term is added to the right hand side (explicit term) of (18) and is of 
the form 

Mo«A— S 2 SyS\ (32) 

where /J. av is of order unity. We found 1.5 < fi av < 2 to be most efficacious. Without the artificial 
viscosity term, the allowed value of A is about one order of magnitude less (e.g., for R = Q = 1, A 
without using AV is about .035). 
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3 COMPACT FOUR STEP FOURTH ORDER RUNGE-KUTTA 
ALGORITHM 

The basic four step explicit Runge-Kutta scheme, as proposed in [2], takes the form: 


u (!) = «(") - 
u( 2 ) = u( n ) — 

„(3) = u («) _ 1 AtR(uW) 
u (*) = u (n) _ AtR(u( 3 )) 

u ( n +l) _ u (4) 


where R is the finite difference (or finite volume) representation of the steady operator. For example, 
in the 2-D second order spatial accuracy case 


*(«) 


M»g» . A 

Ax J Ay 9 


(34) 


where / and g are the flux vectors encountered in section 2. If we want fourth order spatial accuracy, 
then it follows directly that the residual R is modified to read 




and we still retain the compact support. 

Similarly, in the three dimensional case we have 


(35) 
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(36) 


Thus (33), with R given either by (35) for the 2-D case or by (36) for the 3-D case, retains all the 
features of the second order scheme but gains us the fourth order accuracy. In addition one can 
easily verify by simple analysis that for a given cartesian grid and flow conditions the new fourth 
order formulation enhances the stability condition. In the 2-D case we have, using (35) rather than 
(34) 

(At) ( 4 ) ^ (At) 4 th order _ j gg 
(At)( 2 j (At) 2nd order 

In the 3-D case the gain is even more favorable, 


(37) 


(At)4 

(At) 2 


1 . 66 . 


(38) 


Thus the algorithm efficiency gains are two fold. First, for a given acceptable error level the fourth 
order accuracy allows a coarser grid, i.e., fewer node points. Second, not only At is increased due 
to the larger cell size but in addition it gains due to (37) (or (38) in the 3-D case). 
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4 NUMERICAL RESULTS FOR 2-D CASE 


The two-dimensional fourth order compact scheme is applied to a shock reflection problem sketched 
in Figure 1. It shows a 5° shock at Mach 1.95 reflecting from a flat plate. Results are presented 
both without any explicit artificial viscosity (AV) and with a fourth order AV term of the type 
— added to the right hand side. Addition of this AV term reduced the accuracy of the 

compact scheme to third order (note that in 3-D the sixth order AV terms precludes this reduction 
in accuracy). Calculations are also made with a second order implicit Euler AF scheme without and 
with the preceding AV term. Figures 2a and 2b show the results of the compact scheme whereas 
Figures 3a and 3b show the corresponding results from the second order scheme. It is seen that for 
e = 0, both fourth and second order schemes produce very oscillatory results but with e = 0.36, the 
results from the compact scheme (Figure 2b) improve dramatically and, in fact, are much better 
than the corresponding results obtained from the second order scheme (Figure 3b). The shocks 
captured by the compact scheme are sharper and the convergence is also seen to improve. 

Results from a study of grid aspect ratio effect on the compact scheme are also presented here 
in terms of the similarity parameter k of reference 1. Three values of k are considered, namely 1.67, 
1.01, and 0.42. Figure 2b corresponds to k = 1.67 and figures 4 and 5 correspond to k = 1.01 and 

0.42, respectively. In all these cases, e is set equal to 0.36. It is clear from the figures displaying 
the effect of k that the best results are obtained for k near unity. In reference 1, a linear theory 
(e.g., for weak shock) predicts the same results. It is interesting that we find numerically that this 
is also the case in the present nonlinear problem. 

SUMMARY 

1. The steady state solution of the Euler equations of gas dynamics may be achieved to fourth 
order accuracy using a compact grid stencil of 3 x 3 and 3 x 3 X 3 in the 2-D and 3-D cases 
respectively. We presented two examples of such algorithms: one implicit (Euler approximate 
factorizations scheme) and one explicit (Four-stage Runge-Kutta). 

2. Numerical experiments were carried out for the 2-D shock reflection problem, using the im- 
plicit algorithm. Comparisons are made with a corresponding second order scheme. The 
results show that the compact higher order scheme offers marked improvement in both accu- 
racy and convergence rate. 

3. In connection with this work we would like to make the following remarks. It is known that 
the finite difference scheme cannot obtain high order accuracy for conservation-form equations 
computed on a non-uniform grid. This observation coupled with the ease of obtaining fourth 
order compact schemes even in 3-D for the Euler equations revives an old debate: Can one use 
uniform grid and apply conveniently boundary conditions to arbitrary shapes. The potential 
gain in reduced number of computational nodes and enhanced convergence rate appears large 
enough to study this question again. 
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Figure 1: Shock reflection problem 
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